Molecular Dissociation in Hot, Dense Hydrogen 



W. R. Magro 1 , D. M. Ceperley 2 , C. Pierleoni 3 , and B. Bermi 4 
1 Theory Center and Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY 14-853 
2 National Center for Supercomputing Applications and Department of Physics, University of Illinois at Urbana- Champaign, 

Urbana, Illinois 61801 

3 Dipartimento di Fisica, Universita de LAquila 1-67100 LAquila (Italy) and INFM, sezione di Roma I, 00185 Roma 
4 Laboratoire de Physique Theorique des Liquides, Universite P. et M. Curie, 75252 Paris Cedex 05 

(Received February 6, 2008) 

We present a path-integral Monte Carlo study of dissociation in dense hydrogen (1.75 < r s < 2.2, 
with r B the Wigner sphere radius). As temperature is lowered from 10 5 to 5000 K, a molecular 
hydrogen gas forms spontaneously from a neutral system of protons and electrons. At high density, 
r s < 2.0, thermally activated dissociation is accompanied by decreasing pressure, signaling the 
presence of a first order transition and critical point. The decrease in electron kinetic energy during 
dissociation is responsible for the pressure decrease and transition. At lower density the phase 
transition disappears. 
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Despite its simple composition, hydrogen is a complex 
substance with a rich phase diagram. While recent at- 
tention has been directed primarily at its ground state 
properties, much of its behavior at high temperature and 
density, particularly important in astrophysics, remains 
uncertain. Molecular dissociation can occur thermally 
or through compression, as the electrons are forced into 
high energy states by the increasing chemical potential. 
A fundamental question is whether hydrogen crosses a 
phase boundary as the dense molecular fluid transforms 
into a fully ionized plasma. The answer is not simple to 
deduce, since dissociation occurs in a region where both 
thermal and degeneracy effects, as well as many-body 
quantum effects, are important. Further, as in a liquid- 
gas transition, there is no change in symmetry. If such a 
"plasma phase transition" , as it has come to be known, 
indeed exists, it remains unclear whether the dissocia- 
tion and ionization processes occur simultaneously or if 
instead there exists an intermediate atomic-like state. 

Direct observation of the dissociation process is hin- 
dered by the high pressures and temperatures required. 
Recent shock-wave compression measurements, the most 
promising experimental approach, see no evidence for a 
first order phase transition jD. 

One theoretical approach to these questions has been 
to separately model the molecular and metallic phases, 
then equate the Gibbs free energies to locate a transition. 
More advanced models generally include more chemical 
species, such as H2, H, protons, and electrons. A partic- 
ularly sophisticated chemical treatment, the free energy 
model of Saumon and Chabrier , predicts a first order 
phase transition from the molecular fluid to a partially 
ionized atomic gas. The density discontinuity associated 
with this phase transition would sharply alter current 
estimates of the interior mass distribution of the giant 
planets. Chabrier, et al. have shown that a plasma phase 
transition is in fact required to obtain agreement between 
the most sophisticated structural models of Saturn and 
the measured gravitational moments H . 

At such elevated temperatures and densities, it is un- 
clear whether a chemical picture is adequate, since the 
very chemical species in the system are evolving and the 
relevance of the terms 'atom' and 'molecule' is uncer- 
tain. The restricted path-integral Monte Carlo (RPIMC) 
method is unique in its ability to simulate fully interact- 
ing many-fermion quantum systems in thermodynamic 
equilibrium with a minimum of approximations. In par- 
ticular, hydrogen can be modeled within RPIMC as a 
collection of fully interacting electrons and protons. This 
relatively new method has been previously applied only 
to 3 He El, isotopic helium mixtures M, and the dense 
hydrogen plasma ||. We apply it here to investigate 
the nature of molecular dissociation in dense hydrogen, 
thus circumventing the problems associated with chemi- 
cal models. 

We model hydrogen as a neutral mixture of 32 pro- 



tons and 32 unpolarized electrons in a periodically re- 
peated cubic cell and in equilibrium at a temperature, 
T = l/ksf3. Density is specified in terms of the Wigner 
sphere radius, r Sl defined by 47r/3(r s ao) 3 = n , where 
ao is the Bohr radius and n is the average electron den- 
sity. We use the fully interacting, non-relativistic Hamil- 
tonian for this system. The density matrix, p(/3) = e _ ^ w , 
contains complete thermodynamic information about the 
system with observables given as: 

_ TT[Op(f3)} JdR (R\Op((3)\K) 
[ } Tr[p(0)] J dR (R\p((3)\R) ' 1 > 

where R = {r 1; . . . ,rjy} specifies a configuration of the 
N particles. Due to its exponential form, the density 
matrix can be factored as p(f3) = [p(r)] if M = f3/r. 
Eq. (Q) then becomes a path-integral which is well-suited 
for Monte Carlo evaluation using a multi-stage Metropo- 
lis algorithm §. The problem is thus reduced to one of 
evaluating off-diagonal elements of the high-temperature 
density matrix, p(r). 

To develop an expression for p(r), we first split the 
Coulomb potential into one short- and one long-ranged 
term. The high-temperature density matrix, p s (R, R';t), 
in the absence of the long-range potential is written as 
a product of pairwise terms, which are computed with a 
matrix-squaring technique Next, we introduce the 

long-range interaction perturbatively and write 

p(R,R' ] T)=p s (R,R' ] r)e-i[ u ^ +u ^^ (2) 

with the long-range action, C/(i?;r), determined using 
the random phase approximation ^ (details are given in 
Ref. jl(|). Evaluation of the path- integral requires a non- 
zero r, which introduces a systematic 'time-step error,' 
which decreases with decreasing r. In the present work 
r -1 = 10 6 /csK, which gives a reasonable tradeoff be- 
tween computational effort and accuracy. At this value, 
the total energy is accurate to about 0.06 Ry or 5% per 
atom. 

In the path-integral formulation, electrons and pro- 
tons are put on equal footing from a quantum- mechanical 
point of view, since each particle is fully represented 
by a Feynman path. This causes no particular diffi- 
culty, in contrast to other methods, such as local den- 
sity functional calculations. The Fermi temperatures of 
the protons and electrons at the highest density consid- 
ered, r s = 1.75, are 103 K and 190000 K, respectively. 
While the Fermi character of the protons is negligible 
under the conditions of this study, the effects of elec- 
tron exchange are substantial. Without Pauli repulsion, 
a proton-electron mixture is thermodynamically unsta- 
ble flTT| ]. The density matrix must be anti-symmetric 
under exchange of spin- like electrons, but a direct anti- 
symmetrization procedure will be statistically very inef- 
ficient due to cancelation of negative and positive terms. 
This is the well-known sign problem of fermion Monte 
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Carlo techniques. We circumvent this problem by using 
the fixed-node approximation, in which the paths are re- 
stricted to lie within a set of physically motivated trial 
nodes jl2],Q . The procedure becomes exact when the trial 
nodes coincide with the exact fermion nodes. We use the 
nodes of the free-particle density matrix. While these 
nodes become exact only at high temperature, they cap- 
ture more of the physics than one might expect. Their 
merits are discussed in more detail in Ref. In the 
absence of exact results, the errors due to incorrect trial 
nodes can only be estimated from analogous ground state 
calculations Q. These show that for reasonable nodes, 
the energy is largely insensitive to exact nodal positions. 
Typical fixed-node errors at T = are about 0.004 Ry 
per atom, and the error at finite temperature should be 
smaller. 

Dense hydrogen is therefore not only an interesting 
problem, but also an important test case for the method. 
Previously, we successfully applied the same method to a 
dense hydrogen plasma, obtaining good agreement with 
theoretical predictions (|]. We have also tested our pro- 
gram on the isolated hydrogen atom, hydrogen molecule, 
and helium atom. 

In brief, our results for dense hydrogen near disso- 
ciation generally support the findings of Saumon and 
Chabrier from their chemical model Q , though the quan- 
titative details differ. Most significantly, both approaches 
find behaviors consistent with and suggestive of a first or- 
der plasma phase transition. As in the chemical model, 
our molecular gas dissociates first into a partially ion- 
ized atomic-like fluid, then gradually transforms into an 
ionized plasma. 

At the lowest temperatures and densities considered 
(r s = {1.75, 1.86, 2.0, 2.2} , T = 5000 K), a molecular hy- 
drogen gas forms, with the bond somewhat contracted 
from its free space length, 0.742 A. At r s = 2.2 the bond 
length is 0.67 A and further decreases with increasing 
density, reaching 0.65 A at r s = 1.75. The bond con- 
traction is nearly temperature-independent and appar- 
ently results from a stiff effective intermolecular repul- 
sion. This repulsion also leads to an excluded region sur- 
rounding each molecule. From g pp (r), the proton-proton 
pair distribution, at (r s — 2.2, T — 5000 K) we estimate 
the radius of the repulsive core to be 0.6 A, somewhat 
smaller than previous estimates H . 

As the temperature increases, dissociation occurs, as 
is evident from the correlation functions shown in Fig. [l|. 
Dissociation also results from isothermal compression. 
As expected, dn/dT < along the phase boundary, due 
to cooperative thermal and pressure effects on the elec- 
trons. While a first-order dissociation transition would 
proceed isothermally at constant pressure, it occurs in a 
temperature interval at constant volume. At r s — 1.86, 
dissociation occurs for 6000 < T < 8000 K, compared 
with 13000 < T < 15000 K in the chemical model @. 
For Coulomb systems, the pressure is P — (n/3){E + K], 



where E and K are, respectively, the total and kinetic 
energy per atom. Except at r s = 2.2, the pressure de- 
creases during isochoric dissociation, as shown in Fig. ^ 
for r s = 2. In the corresponding isobaric system, this un- 
usual behavior becomes a positive density discontinuity, 
consistent with the negatively sloping phase boundary. 
This is a strong indication of the presence of a first order 
phase transition which terminates in a critical point near 
(r s « 2.2, Tk, 11000 K). 

A possible explanation for the existence of first order 
transition with these behaviors lies in the increasing ki- 
netic energy associated with bond formation. As shown 
in Fig. |^, the electronic kinetic energy normally decreases 
with decreasing temperature. As molecules form, how- 
ever, the kinetic energy increases. The total kinetic en- 
ergy per atom for an isolated atom and H2 molecule are 
also shown for comparison. Since g pe (r) is nearly in- 
variant during molecular formation, the increase in ki- 
netic energy results primarily from angular localization 
as electrons leave spherical atomic-like states in favor of 
molecular bonding states. Bond formation is clearly sig- 
naled by the pairing of spin-unlike electrons, as shown 
in Fig. ^. At low density, the increase in kinetic energy 
during bonding is small relative to the total energy gain, 
so the pressure decreases. As the density increases, how- 
ever, the effects of the inter-molecular repulsion begin 
to dominate, the molecules contract, and the additional 
confinement of the electrons leads to an increasingly high 
kinetic cost of binding. Eventually, AP oc [AE + AK] 
changes sign, and the critical behavior appears. 

To conclusively demonstrate the presence of this phase 
transition, we must show that the above results persist 
in the limits r -> and N -> 00. At t" 1 = 10 6 A: B K, 
three-body and higher order terms in the density matrix 
are significant in dense hydrogen. Our pair-product den- 
sity matrix slightly overestimates the probability for two 
electrons to simultaneously occupy the area between a 
pair of protons. This leads to a slightly shortened bond 
length (by 0.01 A) in H 2 , overestimated kinetic, and cor- 
respondingly underestimated potential energy, due to the 
enhanced electron-proton correlation. Pressures, also af- 
fected, may be overestimated by as much as 15 GPa. 
Finite size effects, quite large in the 32 electron ground 
state, are substantially suppressed at finite temperature, 
due to smearing of the Fermi surface. Calculations of the 
free Fermi gas at T = 5000 K and T — 50000 K indicate 
the finite size errors are smaller than the time-step errors. 

It is interesting to consider the nature of the fluid into 
which the dissociating hydrogen transforms. This phase 
is not yet a plasma, as it retains very strong proton- 
electron correlations. It is perhaps best called a par- 
tially ionized atomic fluid, because the electron binds 
to the proton sufficiently long to prevent other elec- 
trons from approaching. This behavior, apparent from 
r 2 [g pe (i") — 1], gradually disappears as the temperature 
or pressure is raised and the fluid becomes fully ionized. 
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A more precise characterization of these electronic bound 
states and their eventual disappearance can be made with 
the aid of natural orbitals but this requires a sep- 
arate calculation of off-diagonal elements of p(f3). We 
hope to perform this analysis in the future. 

In conclusion, we have simulated dense hydrogen flu- 
ids by assembling a fully interacting collection of protons 
and electrons. The dissociation transformation, relevant 
to interior models of the giant planets, occurs at some- 
what lower temperatures than previous estimates. Wc 
have identified in this system behaviors characteristic of 
a first order phase transition. The phase boundary has 
negative slope, due to the cooperative effects of tempera- 
ture and degeneracy in dissociation. Molecular hydrogen 
dissociates not directly into a plasma, but first into a 
partially ionized atomic fluid. 

The spontaneous formation of molecules in this work, 
in which the only inputs are the Hamiltonian and the 
nodal surface of a free fermion gas, is an important 
and encouraging success for the restricted path-integral 
Monte Carlo method for many-body Fermi systems. All 
aspects of the method, including changing the trial nodes, 
can be improved, so more accurate calculations are un- 
derway. Although new calculations may alter these re- 
sults, we expect the basic findings to remain true. 

A detailed tabulation of the hydrogen equation of state 
at many experimentally inaccessible conditions, needed 
by planetary modelers, will be straightforward to ob- 
tain. Other possible applications include alkali metals, 
the electron-hole liquid, and helium-hydrogen mixtures. 
For the molecular fluid and solid at lower temperatures, 
the errors due to free-particle nodes and time step error 
may become substantial, so more sophisticated nodes and 
high-temperature density matrices will be needed before 
these cases can be studied. 
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length slightly contracted from the free space value. 
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FIG. 2. Computed pressures (open circles) in the present 
work and from the chemical picture (squares) for r s = 2.0. 
Statistical errors are smaller than the symbol size. The 
dashed line is the phase coexistence line transition proposed 
by Saumon and Chabrier Our estimate of the same line 
is given by the thick solid line. The ellipses (• ■ ■) indicate the 
line continues to higher pressure. The triangle at T = is the 
ground state pressure [Q. (Data from the chemical model 
do not have a region of negative dP/dT, since r s — 2.0 lies 
entirely in the supercritical region. The dP/dT < behavior 
does appear at higher densities.) 
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FIG. 4. Pair distribution functions, g ee (r), for spin-unlike 
electrons at r s — 2.0. Temperatures shown are 5000, 7813, 
8927, 10000, 31250, and 125000 K. Spin unlike electrons 
pair to form molecular H2 bonds. Inset: spin- like (dashed) 
and spin-unlike (solid) distribution functions at r 3 = 2.0, 
T = 5000 K. 
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FIG. 3. Electron kinetic energy per atom under isochoric 
cooling from the plasma to molecular gas for three densities: 
r s = 1.86 (filled triangles), r s — 2.0 (crosses), and r s = 2.2 
(open circles). Statistical errors are smaller than the symbol 
size. The kinetic energy increases as electrons pair to form 
bonds. Different isochores cross, since compression suppresses 
molecular formation. The open triangles on the energy axis 
denote the kinetic energy per atom for an isolated hydrogen 
atom and molecule, respectively. 
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